t <- c(rep(0, 100), rep(1, 100))
z <- rep(c(rep(0, 50), rep(1, 50)), 2)
x <- c(rep(0, 45), rep(1, 5), rep(0, 25), rep(1, 25), 
       rep(0, 25), rep(1, 25), rep(0, 25), rep(1, 25))

result.1 <- result.2 <- result.3 <- result.4 <- rep(NA, 1000)
for (i in 1:1000) {
  set.seed(i)
  y.1 <- 1 + x * 1 + z * 1 + t * 1 + x * t * -0.5 + 
    rnorm(200, 0, 0.5)  # Scenario (A.1)
  y.2 <- 1 + x * 1 + z * 1 + t * 1 + x * t * -0.5 + z * t * 2 + 
    rnorm(200, 0, 0.5)  # Scenario (A.2)
  result.1[i] <- lm(y.1 ~ x * t + z)$coefficients[5]  # Model (A.1)
  result.2[i] <- lm(y.1 ~ x * t + z * t)$coefficients[5]  # Model (A.2)
  result.3[i] <- lm(y.2 ~ x * t + z)$coefficients[5]  # Model (A.1)
  result.4[i] <- lm(y.2 ~ x * t + z * t)$coefficients[5]  # Model (A.2)
}

# Table A.3
round(mean(result.1), 3)
round(quantile(result.1, c(0.025, 0.975)), 3)

round(mean(result.2), 3)
round(quantile(result.2, c(0.025, 0.975)), 3)

round(mean(result.3), 3)
round(quantile(result.3, c(0.025, 0.975)), 3)

round(mean(result.4), 3)
round(quantile(result.4, c(0.025, 0.975)), 3)